/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright held by original author
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software; you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation; either version 2 of the License, or (at your
    option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM; if not, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

\*---------------------------------------------------------------------------*/

#include "immersedBoundaryAdjustPhi.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "immersedBoundaryFvPatchFields.H"
#include "immersedBoundaryFvsPatchFields.H"
#include "inletOutletFvPatchFields.H"
#include "fvc.H"

// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //

void Foam::immersedBoundaryAdjustPhi
(
    surfaceScalarField& phi,
    volVectorField& U
)
{
    const fvMesh& mesh = phi.mesh();

    // If the mesh is moving, adjustment needs to be calculated on
    // relative fluxes.  HJ, 13/Feb/2009
    if (mesh.moving())
    {
        fvc::makeRelative(phi, U);
    }

    forAll (phi.boundaryField(), patchI)
    {
        const fvPatchVectorField& Up = U.boundaryField()[patchI];

        if (isA<immersedBoundaryFvPatchVectorField>(Up))
        {
            if (Up.fixesValue())
            {
                // Found immersed boundary path which fixes value.
                // Correction is necessary

                // Cast into immersedBoundary types
                const immersedBoundaryFvPatchVectorField& ibU =
                    refCast<const immersedBoundaryFvPatchVectorField>(Up);

                const immersedBoundaryFvPatch& ibPatch = ibU.ibPatch();

                // Complete immersed boundary triangular surface is present
                // on all processors: no need for parallel reduction

                // Sum the flux through the immersed boundary
                const scalar triFlux = sum(ibU.refValue() & ibPatch.triSf());

                scalarField& phiInternal = phi.internalField();

                scalar fluxIn = 0;
                scalar fluxOut = 0;
                scalar fixedFlux = 0;

                // Get all IB faces
                const labelList& ibFaces = ibPatch.ibFaces();
                const boolList& ibFaceFlips = ibPatch.ibFaceFlips();

                forAll (ibFaces, faceI)
                {
                    const label curFace = ibFaces[faceI];
                    const bool curFlip = ibFaceFlips[faceI];

                    if (mesh.isInternalFace(curFace))
                    {
                        const scalar curFlux = phiInternal[curFace];

                        if (!curFlip)
                        {
                            // Face points out of the live cell
                            if (curFlux >= 0)
                            {
                                // Flux out of the live cell
                                fluxOut += curFlux;
                            }
                            else
                            {
                                // Flux into the live cell
                                fluxIn -= curFlux;
                            }
                        }
                        else
                        {
                            // Face points into the live cell: flip it
                            if (curFlux >= 0)
                            {
                                // Flux into of the live cell
                                fluxIn += curFlux;
                            }
                            else
                            {
                                // Flux out the live cell
                                fluxOut -= curFlux;
                            }
                        }
                    }
                    else
                    {
                        const label patchID =
                            mesh.boundaryMesh().whichPatch(curFace);
                        const label faceID =
                            mesh.boundaryMesh()[patchID].whichFace(curFace);

                        const scalar curFlux =
                            phi.boundaryField()[patchID][faceID];

                        // Note: only coupled patches may carry flux
                        // In order to avoid double summation and
                        // inconsistencies in the correction,
                        // coupled face fluxes will NOT be corrected,
                        // but only accounted for in the summation.
                        if (mesh.boundaryMesh()[patchID].coupled())
                        {
                            // Only do the master side; slave will
                            // be handled on the other side of the couple
                            if (!curFlip)
                            {
                                fixedFlux += curFlux;
                            }
                        }
                    }
                }

                reduce(fluxIn, sumOp<scalar>());
                reduce(fluxOut, sumOp<scalar>());
                reduce(fixedFlux, sumOp<scalar>());

                scalar imbalance = (fluxIn - fluxOut + fixedFlux) - triFlux;

                if (fvMesh::debug)
                {
                    Info<< "triFlux = " << triFlux
                        << " fluxIn = " << fluxIn << " fluxOut = " << fluxOut
                        << " fixedFlux = " << fixedFlux
                        << " imbalance = " << imbalance
                        << endl;
                }

                scalar massCorr = 1.0;

                if (mag(imbalance) > SMALL)
                {
                    // Scaling required: scale to match the smaller of two
                    // fluxes
                    if (fluxIn > fluxOut)
                    {
                        // Scale down incoming flux
                        // Note change of sign: imbalance is negative
                        massCorr = 1 - imbalance/(fluxIn + SMALL);

                        if (fvMesh::debug)
                        {
                            Info<< "Scaling down incoming flux with factor = "
                                << massCorr << endl;
                        }

                        scalar newFluxIn = 0;

                        // Visit all incoming flux faces and re-scale the flux
                        forAll (ibFaces, faceI)
                        {
                            const label curFace = ibFaces[faceI];
                            const bool curFlip = ibFaceFlips[faceI];

                            if (mesh.isInternalFace(curFace))
                            {
                                // Take reference to current flux
                                scalar& curFlux = phiInternal[curFace];

                                if (!curFlip)
                                {
                                    // Face points out of the live cell
                                    if (curFlux < 0)
                                    {
                                        // Flux out of the live cell
                                        curFlux *= massCorr;
                                        newFluxIn -= curFlux;
                                    }
                                }
                                else
                                {
                                    // Face points into the live cell: flip it
                                    if (curFlux >= 0)
                                    {
                                        // Flux out the live cell
                                        curFlux *= massCorr;
                                        newFluxIn += curFlux;
                                    }
                                }
                            }
                        }
                    }
                    else
                    {
                        // Scale down outgoing flux
                        massCorr = 1 + imbalance/(fluxOut + SMALL);

                        if (fvMesh::debug)
                        {
                            Info<< "Scaling down outgoing flux with factor = "
                                << massCorr << endl;
                        }

                        scalar newFluxOut = 0;

                        // Visit all outgoing flux faces and re-scale the flux
                        forAll (ibFaces, faceI)
                        {
                            const label curFace = ibFaces[faceI];
                            const bool curFlip = ibFaceFlips[faceI];

                            if (mesh.isInternalFace(curFace))
                            {
                                // Take reference to current flux
                                scalar& curFlux = phiInternal[curFace];

                                if (!curFlip)
                                {
                                    // Face points out of the live cell
                                    if (curFlux >= 0)
                                    {
                                        // Flux out of the live cell
                                        curFlux *= massCorr;
                                        newFluxOut += curFlux;
                                    }
                                }
                                else
                                {
                                    // Face points into the live cell: flip it
                                    if (curFlux < 0)
                                    {
                                        // Flux out the live cell
                                        curFlux *= massCorr;
                                        newFluxOut -= curFlux;
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }


    // If the mesh is moving, adjustment needs to be calculated on
    // relative fluxes.  Now reverting to absolute fluxes.  HJ, 13/Feb/2009
    if (mesh.moving())
    {
        fvc::makeAbsolute(phi, U);
    }
}


void Foam::immersedBoundaryAdjustPhi
(
    surfaceScalarField& phi,
    volVectorField& U,
    const surfaceScalarField& phiU
)
{
    const fvMesh& mesh = phi.mesh();

    // If the mesh is moving, adjustment needs to be calculated on
    // relative fluxes.  HJ, 13/Feb/2009
    if (mesh.moving())
    {
        fvc::makeRelative(phi, U);
    }

    forAll (phi.boundaryField(), patchI)
    {
        const fvPatchVectorField& Up = U.boundaryField()[patchI];

        if (isA<immersedBoundaryFvPatchVectorField>(Up))
        {
            if (Up.fixesValue())
            {
                // Found immersed boundary path which fixes value.
                // Correction is necessary

                // Cast into immersedBoundary types
                const immersedBoundaryFvPatchVectorField& ibU =
                    refCast<const immersedBoundaryFvPatchVectorField>(Up);

                const immersedBoundaryFvPatch& ibPatch = ibU.ibPatch();

                // Complete immersed boundary triangular surface is present
                // on all processors: no need for parallel reduction

                // Sum the flux through the immersed boundary
                const scalar triFlux = sum(ibU.refValue() & ibPatch.triSf());

                scalarField& phiInternal = phi.internalField();

                scalar fluxIn = 0;
                scalar fluxOut = 0;
                scalar fixedFlux = 0;

                // Get all IB faces
                const labelList& ibFaces = ibPatch.ibFaces();
                const boolList& ibFaceFlips = ibPatch.ibFaceFlips();

                forAll (ibFaces, faceI)
                {
                    const label curFace = ibFaces[faceI];
                    const bool curFlip = ibFaceFlips[faceI];

                    if (mesh.isInternalFace(curFace))
                    {
                        phiInternal[curFace] = phiU.internalField()[curFace];

                        const scalar curFlux = phiInternal[curFace];

                        if (!curFlip)
                        {
                            // Face points out of the live cell
                            if (curFlux >= 0)
                            {
                                // Flux out of the live cell
                                fluxOut += curFlux;
                            }
                            else
                            {
                                // Flux into the live cell
                                fluxIn -= curFlux;
                            }
                        }
                        else
                        {
                            // Face points into the live cell: flip it
                            if (curFlux >= 0)
                            {
                                // Flux into of the live cell
                                fluxIn += curFlux;
                            }
                            else
                            {
                                // Flux out the live cell
                                fluxOut -= curFlux;
                            }
                        }
                    }
                    else
                    {
                        const label patchID =
                            mesh.boundaryMesh().whichPatch(curFace);
                        const label faceID =
                            mesh.boundaryMesh()[patchID].whichFace(curFace);

                        const scalar curFlux =
                            phi.boundaryField()[patchID][faceID];

                        // Note: only coupled patches may carry flux
                        // In order to avoid double summation and
                        // inconsistencies in the correction,
                        // coupled face fluxes will NOT be corrected,
                        // but only accounted for in the summation.
                        if (mesh.boundaryMesh()[patchID].coupled())
                        {
                            // Only do the master side; slave will
                            // be handled on the other side of the couple
                            if (!curFlip)
                            {
                                fixedFlux += curFlux;
                            }
                        }
                    }
                }

                reduce(fluxIn, sumOp<scalar>());
                reduce(fluxOut, sumOp<scalar>());
                reduce(fixedFlux, sumOp<scalar>());

                scalar imbalance = (fluxIn - fluxOut + fixedFlux) - triFlux;

                if (fvMesh::debug)
                {
                    Info<< "triFlux = " << triFlux
                        << " fluxIn = " << fluxIn << " fluxOut = " << fluxOut
                        << " fixedFlux = " << fixedFlux
                        << " imbalance = " << imbalance
                        << endl;
                }

                scalar massCorr = 1.0;

                if (mag(imbalance) > SMALL)
                {
                    // Scaling required: scale to match the smaller of two
                    // fluxes
                    if (fluxIn > fluxOut)
                    {
                        // Scale down incoming flux
                        // Note change of sign: imbalance is negative
                        massCorr = 1 - imbalance/(fluxIn + SMALL);

                        if (fvMesh::debug)
                        {
                            Info<< "Scaling down incoming flux with factor = "
                                << massCorr << endl;
                        }

                        scalar newFluxIn = 0;

                        // Visit all incoming flux faces and re-scale the flux
                        forAll (ibFaces, faceI)
                        {
                            const label curFace = ibFaces[faceI];
                            const bool curFlip = ibFaceFlips[faceI];

                            if (mesh.isInternalFace(curFace))
                            {
                                // Take reference to current flux
                                scalar& curFlux = phiInternal[curFace];

                                if (!curFlip)
                                {
                                    // Face points out of the live cell
                                    if (curFlux < 0)
                                    {
                                        // Flux out of the live cell
                                        curFlux *= massCorr;
                                        newFluxIn -= curFlux;
                                    }
                                }
                                else
                                {
                                    // Face points into the live cell: flip it
                                    if (curFlux >= 0)
                                    {
                                        // Flux out the live cell
                                        curFlux *= massCorr;
                                        newFluxIn += curFlux;
                                    }
                                }
                            }
                        }
                    }
                    else
                    {
                        // Scale down outgoing flux
                        massCorr = 1 + imbalance/(fluxOut + SMALL);

                        if (fvMesh::debug)
                        {
                            Info<< "Scaling down outgoing flux with factor = "
                                << massCorr << endl;
                        }

                        scalar newFluxOut = 0;

                        // Visit all outgoing flux faces and re-scale the flux
                        forAll (ibFaces, faceI)
                        {
                            const label curFace = ibFaces[faceI];
                            const bool curFlip = ibFaceFlips[faceI];

                            if (mesh.isInternalFace(curFace))
                            {
                                // Take reference to current flux
                                scalar& curFlux = phiInternal[curFace];

                                if (!curFlip)
                                {
                                    // Face points out of the live cell
                                    if (curFlux >= 0)
                                    {
                                        // Flux out of the live cell
                                        curFlux *= massCorr;
                                        newFluxOut += curFlux;
                                    }
                                }
                                else
                                {
                                    // Face points into the live cell: flip it
                                    if (curFlux < 0)
                                    {
                                        // Flux out the live cell
                                        curFlux *= massCorr;
                                        newFluxOut -= curFlux;
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }


    // If the mesh is moving, adjustment needs to be calculated on
    // relative fluxes.  Now reverting to absolute fluxes.  HJ, 13/Feb/2009
    if (mesh.moving())
    {
        fvc::makeAbsolute(phi, U);
    }
}

// ************************************************************************* //
